## vector with package names
x <- c("pbapply", "viridis", "knitr", "kableExtra", "ggplot2", "ape", "PhenotypeSpace", "lmerTest", "brms")
aa <- lapply(x, function(y) {
# check if installed, if not then install
if (!y %in% installed.packages()[,"Package"])
install.packages(y)
# load package
try(require(y, character.only = T), silent = T)
})
knitr::opts_knit$set(root.dir = normalizePath(".."))
knitr::opts_chunk$set(dpi = 38, fig.width = 18, fig.height = 10, echo = TRUE)
options(knitr.kable.NA = '')
source('~/Dropbox/Projects/lbh_cultural_evolution_revBayes/source/custom_treespace.R')
theme_set(theme_classic(base_size = 22))
rb.trees <- readRDS("./output/revbayes_output_in_single_R_object.RDS")
rb.trees <- rb.trees[grep("run_", names(rb.trees), invert = TRUE, value = TRUE)]
trees_diags <- data.frame(model = names(rb.trees))
# get lek name
trees_diags$lek <- sapply(strsplit(trees_diags$model, "_", fixed = TRUE), "[", 1)
# substitution model
trees_diags$subs <- sapply(strsplit(trees_diags$model, "_", fixed = TRUE), "[", 2)
# period
trees_diags$period <- sapply(strsplit(trees_diags$model, "_", fixed = TRUE), "[", 3)
# and which fossils were used
trees_diags$fossils <- sapply(strsplit(trees_diags$model, "_", fixed = TRUE), "[", 4)
trees_diags$align <- sapply(strsplit(trees_diags$model, "_", fixed = TRUE), "[", 5)
trees_diags$n_trees <- sapply(rb.trees, function(x) length(x$trees))
# Number of models by iterations and lek
kbl <- knitr::kable(as.matrix(table(trees_diags$n_trees, trees_diags$lek)), caption = "Number of models with a specific number of trees by lek")
kableExtra::kable_styling(kbl)
# combine data sets
topo_size_all$fossils <- "all"
topo_size_early$fossils <- "early"
topo_size <- rbind(topo_size_all, topo_size_early)
iterations <- 5000
priors <- c(set_prior("normal(0, 1.5)", class = "Intercept"), set_prior("normal(0, 1.5)", class = "b"))
size_mod <- brm(size ~ alignment + data.set + models + fossils, data = topo_size, prior = priors, iter = iterations, chains = 3, silent = 2, seed = 5)
priors <- c(set_prior("normal(0, 1.5)", class = "b"))
# without intercept
size_mod_no_intercept <- brm(size ~ alignment + data.set + models + fossils -1, data = topo_size, prior = priors, iter = iterations, chains = 3, silent = 2, seed = 5)
saveRDS(list(size_mod = size_mod, size_mod_no_intercept = size_mod_no_intercept), "./data/processed/regression_results_topological_size.RDS")
Intercept-included model:
size_mods <- readRDS("./data/processed/regression_results_topological_size.RDS")
size_mods$size_mod
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: size ~ alignment + data.set + models + fossils
## Data: topo_size (Number of observations: 96)
## Draws: 3 chains, each with iter = 5000; warmup = 2500; thin = 1;
## total post-warmup draws = 7500
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.32 0.08 0.16 0.49 1.00 11144 5436
## alignmentoptimal -0.00 0.09 -0.18 0.17 1.00 9896 5775
## alignmentprank -0.12 0.09 -0.30 0.05 1.00 9964 6014
## data.setold 0.35 0.07 0.21 0.49 1.00 12043 5663
## modelsUexp 0.05 0.07 -0.08 0.19 1.00 10779 5259
## fossilsearly -0.00 0.07 -0.14 0.14 1.00 12679 5083
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.34 0.03 0.30 0.40 1.00 10381 5511
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Contrasts for alignment strategies:
aligment_contrasts <- c(prnk_vs_agnostic = "alignmentprank - alignmentall.equal = 0", prnk_vs_optimal = "alignmentprank - alignmentoptimal = 0", optimal_vs_agnostic = "alignmentoptimal - alignmentall.equal = 0")
hypothesis(size_mods$size_mod_no_intercept, aligment_contrasts)
## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
## 1 prnk_vs_agnostic -0.12 0.09 -0.29 0.05 NA NA
## 2 prnk_vs_optimal -0.12 0.09 -0.29 0.05 NA NA
## 3 optimal_vs_agnostic 0.00 0.09 -0.17 0.16 NA NA
## Star
## 1
## 2
## 3
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
topo_similarity_all$fossils <- "all"
topo_similarity_early$fossils <- "early"
topo_similarity <- rbind(topo_similarity_all, topo_similarity_early)
similarity_mod <- brm(overlap ~ alignment + data.set + models + fossils, data = topo_similarity, prior = priors, iter = iterations, chains = 3, silent = 2, seed = 5)
# without intercept
similarity_mod_no_intercept <- brm(overlap ~ alignment + data.set + models + fossils -1, data = topo_similarity, iter = iterations, chains = 3, silent = 2, seed = 5)
saveRDS(list(similarity_mod = similarity_mod, similarity_mod_no_intercept = similarity_mod_no_intercept), "./data/processed/regression_results_topological_similiarity.RDS")
Intercept-included model:
similarity_mods <- readRDS("./data/processed/regression_results_topological_similiarity.RDS")
similarity_mods$similarity_mod
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: overlap ~ alignment + data.set + models + fossils
## Data: topo_similarity (Number of observations: 912)
## Draws: 3 chains, each with iter = 5000; warmup = 2500; thin = 1;
## total post-warmup draws = 7500
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.15 0.02 0.12 0.18 1.00 9958 5540
## alignmentoptimal -0.03 0.02 -0.07 0.00 1.00 9358 5836
## alignmentprank -0.09 0.02 -0.13 -0.06 1.00 9072 6059
## data.setold -0.11 0.01 -0.14 -0.08 1.00 11548 5900
## modelsUexp -0.00 0.01 -0.03 0.03 1.00 10243 5235
## fossilsearly 0.00 0.01 -0.03 0.03 1.00 11138 5780
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.22 0.01 0.21 0.23 1.00 11619 5798
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Contrasts for alignment strategies:
hypothesis(similarity_mods$similarity_mod_no_intercept, aligment_contrasts)
## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob
## 1 prnk_vs_agnostic -0.09 0.02 -0.13 -0.06 NA NA
## 2 prnk_vs_optimal -0.06 0.02 -0.10 -0.03 NA NA
## 3 optimal_vs_agnostic -0.03 0.02 -0.07 0.00 NA NA
## Star
## 1 *
## 2 *
## 3
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
R session information
## R version 4.1.1 (2021-08-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
##
## locale:
## [1] LC_CTYPE=es_ES.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=es_CR.UTF-8 LC_COLLATE=es_ES.UTF-8
## [5] LC_MONETARY=es_CR.UTF-8 LC_MESSAGES=es_ES.UTF-8
## [7] LC_PAPER=es_CR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] brms_2.16.3 Rcpp_1.0.8 lmerTest_3.1-3
## [4] lme4_1.1-27.1 Matrix_1.3-4 PhenotypeSpace_0.1.0
## [7] warbleR_1.1.27 NatureSounds_1.0.4 seewave_2.1.8
## [10] tuneR_1.3.3.1 ape_5.5 ggplot2_3.3.5
## [13] kableExtra_1.3.4 knitr_1.37 viridis_0.6.2
## [16] viridisLite_0.4.0 pbapply_1.5-0
##
## loaded via a namespace (and not attached):
## [1] backports_1.2.1 fastmatch_1.1-3 systemfonts_1.0.2
## [4] plyr_1.8.6 igraph_1.2.6 sp_1.4-6
## [7] splines_4.1.1 crosstalk_1.1.1 inline_0.3.19
## [10] rstantools_2.1.1 digest_0.6.29 htmltools_0.5.2
## [13] rsconnect_0.8.24 fansi_1.0.2 magrittr_2.0.2
## [16] checkmate_2.0.0 tensor_1.5 cluster_2.1.2
## [19] RcppParallel_5.1.4 matrixStats_0.60.1 xts_0.12.1
## [22] svglite_2.0.0 spatstat.sparse_2.1-0 CircStats_0.2-6
## [25] prettyunits_1.1.1 colorspace_2.0-2 signal_0.7-7
## [28] rvest_1.0.1 xfun_0.29 dplyr_1.0.7
## [31] callr_3.7.0 crayon_1.4.2 RCurl_1.98-1.5
## [34] jsonlite_1.7.2 spatstat.data_2.1-2 phangorn_2.8.1
## [37] zoo_1.8-9 glue_1.6.1 polyclip_1.10-0
## [40] gtable_0.3.0 webshot_0.5.2 distributional_0.2.2
## [43] pkgbuild_1.2.0 rstan_2.21.3 abind_1.4-5
## [46] scales_1.1.1 mvtnorm_1.1-2 DBI_1.1.1
## [49] miniUI_0.1.1.1 dtw_1.22-3 xtable_1.8-4
## [52] diffobj_0.3.4 spatstat.core_2.3-2 proxy_0.4-26
## [55] StanHeaders_2.21.0-7 stats4_4.1.1 DT_0.19
## [58] htmlwidgets_1.5.3 httr_1.4.2 threejs_0.3.3
## [61] posterior_1.0.1 ellipsis_0.3.2 pkgconfig_2.0.3
## [64] loo_2.4.1 farver_2.1.0 sass_0.4.0
## [67] deldir_1.0-6 utf8_1.2.2 labeling_0.4.2
## [70] tidyselect_1.1.1 rlang_1.0.0 reshape2_1.4.4
## [73] later_1.3.0 munsell_0.5.0 tools_4.1.1
## [76] cli_3.1.0 generics_0.1.0 ade4_1.7-18
## [79] adehabitatHR_0.4.19 ggridges_0.5.3 evaluate_0.14
## [82] stringr_1.4.0 fastmap_1.1.0 yaml_2.2.2
## [85] goftest_1.2-3 processx_3.5.2 purrr_0.3.4
## [88] nlme_3.1-152 projpred_2.0.2 mime_0.11
## [91] xml2_1.3.2 compiler_4.1.1 bayesplot_1.8.1
## [94] shinythemes_1.2.0 rstudioapi_0.13 gamm4_0.2-6
## [97] spatstat.utils_2.3-0 tibble_3.1.6 bslib_0.2.5.1
## [100] stringi_1.7.6 highr_0.9 ps_1.6.0
## [103] Brobdingnag_1.2-6 rgeos_0.5-9 lattice_0.20-44
## [106] nloptr_1.2.2.2 markdown_1.1 shinyjs_2.0.0
## [109] vegan_2.5-7 fftw_1.0-6.1 permute_0.9-5
## [112] tensorA_0.36.2 vctrs_0.3.8 pillar_1.6.4
## [115] lifecycle_1.0.1 spatstat.geom_2.3-1 jquerylib_0.1.4
## [118] bridgesampling_1.1-2 bitops_1.0-7 raster_3.5-11
## [121] httpuv_1.6.2 R6_2.5.1 promises_1.2.0.1
## [124] gridExtra_2.3 codetools_0.2-18 boot_1.3-28
## [127] colourpicker_1.1.0 MASS_7.3-54 gtools_3.9.2
## [130] assertthat_0.2.1 rjson_0.2.21 withr_2.4.3
## [133] shinystan_2.5.0 adehabitatMA_0.3.14 mgcv_1.8-36
## [136] parallel_4.1.1 terra_1.5-12 quadprog_1.5-8
## [139] grid_4.1.1 rpart_4.1-15 adehabitatLT_0.3.25
## [142] coda_0.19-4 minqa_1.2.4 rmarkdown_2.10
## [145] numDeriv_2016.8-1.1 shiny_1.6.0 base64enc_0.1-3
## [148] dygraphs_1.1.1.6